link <- "/data/Dagobah/greengrass/grassland_ger/01_CTM_GL_mask_17-19/GL_mask_2017-2019.tif"
mask_GL <- terra::rast(link)
plot(mask_GL)# Full mosaic
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt.ovr"
# base case
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.vrt.ovr"
# Single tile base case
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif"
tsi_base_raw <- rast(link)
# Single tile case day_16_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
tsi_day_16_sigma_8_16_32_raw <- rast(link)
# Single tile case day_8_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_8_sigma_8_16_32/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
tsi_day_8_sigma_8_16_32_raw <- rast(link)
# Single tile case day_16_sigma_4_8_16
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_4_8_16/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
tsi_day_16_sigma_4_8_16_raw <- rast(link)
tsi_list <- list(tsi_base_raw, tsi_day_16_sigma_8_16_32_raw, tsi_day_16_sigma_4_8_16_raw, tsi_day_8_sigma_8_16_32_raw)
tsi_obs_length <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_obs_length[i] <<- length(names(tsi_list[[i]]))
)# Set a "plan" for how the code should run.
plan(multisession, workers = 4)
# This does run in parallel!
future_map(seq_along(tsi_list), function(i)
terra::plot(tsi_list[[i]][[10:10]] )
)
lapply(seq_along(tsi_list),
function(i) plot(tsi_list[[i]][[10:18]])
)lapply(seq_along(tsi_list),
function(i)
tsi_list[[i]] <<- setMinMax(tsi_list[[i]])
)## [[1]]
## class : SpatRaster
## dimensions : 1000, 1000, 765 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4076026, 4106026, 3134920, 3164920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif
## names : 19840~LND05, 19840~LND05, 19840~LND05, 19841~LND05, 19850~LND05, 19850~LND05, ...
## min values : -6318, 1235, -4413, 2248, -1832, -2821, ...
## max values : 8778, 8231, 8899, 7515, 7072, 8358, ...
##
## [[2]]
## class : SpatRaster
## dimensions : 1000, 1000, 866 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4076026, 4106026, 3134920, 3164920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif
## names : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ...
## min values : NaN, NaN, NaN, NaN, -6318, -6318, ...
## max values : NaN, NaN, NaN, NaN, 8778, 8778, ...
##
## [[3]]
## class : SpatRaster
## dimensions : 1000, 1000, 866 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4076026, 4106026, 3134920, 3164920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif
## names : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ...
## min values : NaN, NaN, NaN, NaN, NaN, NaN, ...
## max values : NaN, NaN, NaN, NaN, NaN, NaN, ...
##
## [[4]]
## class : SpatRaster
## dimensions : 1000, 1000, 866 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 4076026, 4106026, 3134920, 3164920 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif
## names : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ...
## min values : NaN, NaN, NaN, NaN, -6318, -6318, ...
## max values : NaN, NaN, NaN, NaN, 8778, 8778, ...
# create filter out all layers that dont contain any values
tsi_na_index <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_na_index[i] <<- minmax(tsi_list[[i]]) |>
as.data.frame() |>
slice(1) |>
pivot_longer(everything()) |>
mutate(id = row_number()) |>
na.omit()
)
# convert into index vector
lapply(seq_along(tsi_list),
function(i)
tsi_na_index[[i]] <<- as.vector(tsi_na_index[[i]])
)
# subset tsi object with only layers containing values
lapply(seq_along(tsi_list),
function(i)
tsi_list[[i]] <<- tsi_list[[i]][[ tsi_na_index[[i]] ]]
)
# filter out reflectance values below and above 0-100 % range
lapply(seq_along(tsi_list),
function(i)
tsi_list[[i]] <- tsi_list[[i]] |>
filter(across(1:length(names(tsi_list[[i]])), ~ . > 0),
across(1:length(names(tsi_list[[i]])), ~ . < 10000))
)lapply(seq_along(tsi_list),
function(i)
tsi_obs_length[[i]] - length(names(tsi_list[[i]]))
)
lapply(seq_along(tsi_list),
function(i)
print(paste( round(( tsi_obs_length[i] / length(names(tsi_list[[i]]))), digits=4), "% of layers consisted completely of NA values (n=", ( tsi_obs_length[i] - length(names(tsi_list[[i]]))), "out of", tsi_obs_length[i], ")"))
)link_a <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0047/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
link_b <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0052/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
a <- rast(link_a)
b <- rast(link_b)
ab <- merge(a,b)
plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
range=c(0,10000)
)
object.size((ab))
writeRaster(ab, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif",
# overwrite=TRUE,
datatype = "INT4S",
filetype='GTiff')
link_ab <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif"
ab <- rast(link_ab)
plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
range=c(0,10000)
)#rasterVis::levelplot(tsi[[c(24)]], par.settings = viridisTheme)
#rasterVis::levelplot(tsi[[c(1,4,24,25:35)]], par.settings = viridisTheme)
lapply(seq_along(tsi_list),
function(i)
plot(tsi_list[[i]][[c(1:10)]], legend=F, col = viridis(n=100,option="D"),
range=c(0,10000)
)
)## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
for (i in c(1:100
#length(names(tsi))
)) {
# Step 1: Call the pdf command to start the plot
png(file = paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/plot_", i, ".png"), width = 265, height = 125, units='mm', res = 100) #, # The directory you want to save the file in
# width = 12, # The width of the plot in inches
#height = 12) # The height of the plot in inches
plot(tsi[[i]], legend=T, col = viridis(n=100,option="D"),
range=c(0,10000),
plg=list( title = paste(as.Date(names(tsi[[i]]), format="%Y%m%d")), title.cex=1.25)
)
# Step 3: Run dev.off() to create the file
dev.off()
}knitr::include_graphics("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/ts_16_8_16_32.gif")tsi_df <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_df[[i]] <<- terra::as.data.frame(tsi_list[[i]][[c(1:100)]], xy = TRUE, na.rm=F) |>
# filter out rows that only contain na
filter_at(vars(-x,-y), any_vars(! is.na(.))) |>
# random subset for data vis
slice_sample(n=10)
)
# create long version of df with unique and formatted location and date cols
tsi_df_long <- list()
lapply(seq_along(tsi_list),
function(i)
tsi_df_long[[i]] <<- tsi_df[[i]] |>
pivot_longer(cols = c(3:ncol(tsi_df[[i]])), names_to = "timestamp") |>
mutate(xy= paste0(x,"_",y)) |>
mutate(date = as.Date(timestamp, format="%Y%m%d"),
year = lubridate::year(date))
)lapply(seq_along(tsi_list),
function(i)
ggplot(tsi_df_long[[i]], aes(date, value, group=xy, color=xy)) +
geom_line() +
geom_point(size=0.5) +
theme_classic() +
scale_colour_viridis_d(option = "D") +
scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position="none") +
facet_wrap(~year, scales = "free_x")
)## [[1]]
## Warning: Removed 42 row(s) containing missing values (geom_path).
## Warning: Removed 434 rows containing missing values (geom_point).
##
## [[2]]
## Warning: Removed 21 row(s) containing missing values (geom_path).
## Warning: Removed 177 rows containing missing values (geom_point).
##
## [[3]]
## Warning: Removed 18 row(s) containing missing values (geom_path).
## Warning: Removed 259 rows containing missing values (geom_point).
##
## [[4]]
## Warning: Removed 42 row(s) containing missing values (geom_path).
## Warning: Removed 247 rows containing missing values (geom_point).
ggplot(tsi_df_long, aes(date, value, group=xy, color=xy)) +
geom_line() +
theme_classic() +
scale_colour_viridis_d(option = "D") +
scale_x_date(date_breaks = "1 month", date_labels = "%m-%y")+
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position="none")
ggplot(tsi_df_long, aes(date, value)) +
geom_smooth() +
theme_classic() +
scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
theme(axis.text.x = element_text(angle = 45 , hjust=1))
ggplot(tsi_df_long, aes(date, value, group = timestamp)) +
geom_boxplot() +
theme_classic() +
scale_x_date( date_labels = "%d-%m-%y")+
theme(axis.text.x = element_text(angle = 45, hjust=1)) ggplot() +
geom_spatraster(data = tsi[[c(4:8)]]) +
facet_wrap(~lyr, ncol = 2) +
scale_fill_whitebox_c(
palette = "viridi",
# labels = scales::label_number(suffix = "ยบ")
) +
labs(fill = "NDVI") +
theme_classic()